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One of the many remarkable properties of graphene is that in the low energy limit the dynamics 
of its electrons can be effectively described by the massless Dirac equation. This has prompted 
investigations of graphene based on the lattice simulation of a system of 2-dimensional fermions 
on a square staggered lattice. We demonstrate here how to construct the path integral for graphene 
working directly on the graphene hexagonal lattice. For the nearest neighbor tight binding model 
interacting with a long range Coulomb interaction between the electrons, this leads to the hybrid 
Monte Carlo algorithm with no sign problem. The only approximation is the discretization of the 
Euclidean time. So as we extrapolate to the time continuum limit, the exact tight binding solution 
maybe found numerically to arbitrary precession on a finite hexagonal lattice. The potential for 
this approach is tested on a single hexagonal cell. 
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1. Introduction 

During the last few years graphene, a single layer of carbon atoms forming a hexagonal lattice, 
has burst on the scientific scene as a remarkable material with intriguing theoretical properties and 
potentially astounding technological applications [1]. Its electronic structure is well approximated 
by the tight binding Hamiltonian H which describes the quantum mechanical motion of the elec- 
trons, one per atom at half-filling, that are not part of the cloud of electrons responsible for the 
covalent binding of the carbon atoms. The simplest model Hamiltonian, H = H2 + He, consists 
of two terms, a quadratic Hamiltonian H2 which describes the hopping of the electrons between 
neighboring atoms, and a Hamiltonian He which describes the Coulomb repulsion of the electrons. 
In the the limit, where the strength of the hopping term is taken to be much larger than strength of 
the Coulomb interaction, the quadratic Hamiltonian H2 gives origin to a dispersion formula which 
for low momenta is analogous to the dispersion formula for relativistic fermions in two dimen- 
sions [2]. This has prompted some researchers to adapt to graphene lattice gauge theory techniques 
which have been profitably used for the study of Quantum Chromodynamics and other particle sys- 
tems [3]. In the work of [4, 5] one approximates the tight-binding Hamiltonian of graphene with 
a Dirac Hamiltonian, incorporates the Coulomb interaction through the introduction of a suitable 
electromagnetic field, and finally discretized the resulting continuum quantum field theory on a 
hypercubic space-time lattice. The hybrid Monte Carlo method [6], widely used in lattice gauge 
theory to simulate fermions interacting with quantum gauge fields, can then be used to investigate 
the graphene system. 

The approach outlined above has led to very interesting and valuable results [4, 5], yet it has 
several limitations. First the use of staggered fermion on a cubic lattice looses the direct connection 
to special symmetries of hexagonal lattice and introduces symmetry breaking artifacts, know in 
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lattice QCD as taste splitting [7]. Second a 4-d lattice is required to simulate the coulomb potential 
which is not only computational expensive but distorted at the scale of the physical spatial lattice 
spacing. One would think that, since the starting point is a system already defined on a lattice, 
it should be possible to apply the hybrid Monte Carlo technique directly to the graphene lattice. 
The clear advantage of working on the hexagonal graphene lattice is the direct connection to the 
experimentally determined physical lattice constants of the tight-binding model, which represents 
an accurate description of the experimental system. Here we wish to illustrate how this can be 
done. 



2. The graphene lattice 



Figure 1: The underlying triangular lattice. Figure 2: The hexagonal graphene lattice. 

It is useful to review the geometrical properties of graphene. Graphene is a system of interact- 
ing electrons located at the vertices of a hexagonal lattice. It is convenient to think of the graphene 
lattice as consisting of two triangular sublattices, which we denote by A and B, which together with 
the centers of the hexagons (sublattice C) form a finer, underlying triangular lattice. In Figures 1 
and 2, we illustrate the full underlying triangular lattice and the hexagonal lattice which one ob- 
tains by eliminating the sites of one of the three sublattices. A and B are Bravais lattices so given 
any two primitive vectors a\,a2 their sites ( at = x\U\ + X2«2 and 7b = x[ai +x' 2 a2 respectively) 
are conveniently enumerated by two integers (jq,^). 

One may impose periodic boundary conditions (PBC) by identifying sites that differ by a 
translation by L,- along the primitive vectors: x,- ~ x,- +I4. In Figure 3, we show an example of PBC 
for L\ = L2 = 4 with green and purple color for the independent lattice sites on each sublattice. 
In this representation, each sublattice A and B is mapped into a rectangle with periodic boundary 
condition very easily coded by two dimensional arrays. The pattern of these sites does not look 
reminiscent of the hexagonal structure of graphene. However, it is equivalent. If one rearranges the 
sites by replacing some of the vertices with their periodic images and one adds to these vertices a 
few of their periodic images one obtains a pattern of sites which clearly shows the hexagonal struc- 
ture of graphene with 3 fold rotational symmetry as illustrated in Figure 4. Of course, one could 
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Figure 3: Periodic lattice with 2L\Lj_ sites viewed 
as parallelogram subdomain in the plane illustrated 
for period L, = 4 along two primitive vectors: d\ ~ 
\^3x/2 +y/2, &2 ~ y, defined relative to the Carte- 
sian unit vectors, (x,y). 



Figure 4: The sublattice with period L = 4 from 
Figure 3 is mapped into periodic hexagonal do- 
main by moving four vertices to equivalent sites 
combined with ten image or ghost sites shown in 
lighter color to complete the pattern. 



also introduce different periods L\,L 2 along two directions of the fundamental cell and consider 
other types of boundary conditions. For example, by taking L 2 3> L\ with PBC along the shortest 
dimension and open boundary conditions direction along the other one can simulate a nanotube. 



3. The dynamics 

We introduce fermionic annihilation and creation operators a x , s ,a[ s for the electrons on the 
two sublattices, where x is a site index for both the A and B sublattices and s = ±1 labels the spin 
of the electrons. The Hamiltonian H = H2+ He consists of two terms. The quadratic kinetic term 
is 

H 2 = £ -jk(4/V +h.c), (3.1) 

(x,y) ,s 

where the sum runs over all pairs (x,y) of nearest neighbor sites (coupling the A and B sublattices) 
and the two values of the spin. The Coulomb interaction term is 
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He = IVV^4 y = - £ q x 
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q.x—qx 

x r 



(3.2) 



where q x = a\ ^a x ,\ +a\ _ Y a x -i — 1 is a local charge operator and V is the interaction potential with 
physical co-ordinates r x ,r y for each Carbon atom. Note that we have introduced a -1 in the charge 
operator q x to account for the background charge of the carbon ion: this ensures that the system 
is neutral at half filling, and it will play an important role for our functional integral formulation. 
Insofar as the matrix V in Eq. 3.2 is concerned, it could be the actual 3d Coulomb potential, but 
could be any other interaction potential. The only constraint is that the matrix V x>y must be positive 
definite. 
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The phenomenological constants are a nearest neighbor lattice spacing of a ~ 1.42, K ~ 
l.SMev and e 2 is effective parameter for e 2 /hv ~ 300 xe 2 /hc so that unlike free electrons effective 
charge is very strong and the electrodynamics essential static, free of magnetic effects. The single 
site Coulomb [8] radius, ro ~ 0.5a, is essential for stability of the vacuum configuration. We also 
note that in Eq. 3.1 we have neglected the smaller next-to-nearest neighbor hopping within each 
sublattice , estimated to be K*' ~ 0.03k". This would introduce a (probably manageable) complex 
phase in the path integral. 

Finally, we observe that the Hamiltonian of Eqs. 3.1-3.2 commutes with the angular momen- 
tum generators, 

J± = a l,s a ± a x,s'i 

J 3 = al s of'a x . s >/2 , (3.3) 

which act on the spin of the electron and play the role in the tight-binding model of an internal 
or "isospin" symmetry. Additional symmetries of the system are the overall momentum conserva- 
tion, parity under reflection with respect to the crystal axes and under interchange of the A and B 
sublattices, and the separate conservations of the number of electrons with spin up and spin down. 
In order to explore the properties of the system one would like to calculate expectation values 

(Mti) 02(h)-..) =Z- { TvT[0 l {h)0 2 {t 2 )e-P H ] , (3.4) 

where & stands for a generic observable, j3 can be interpreted as an extent in Euclidean time, T[...] 
stands for time ordering of the operators inside the square bracket with respect to the Euclidean 
evolution implemented by exp(— fiH), and Z = Tr exp(— f)H) is the partition function. 

4. Path integral formulation 

Our goal is to provide an equivalent path integral formulation of Eq. 3.4 conducive to calcu- 
lation by numerical simulation, following a rather standard procedure to convert from the Hamil- 
tonian into a Lagrangian. We will first express the expectation values and the partition function 
in terms of an integral over anticommuting fermionic fields, i.e. elements of a Grassmann algebra. 
(See for example [9].) It is important that the Hamiltonian be normal ordered. This is true of H 2 
and of the non-local terms in He- However the normal ordered : He : differs from He by a diagonal 
quadratic term which can be added to H 2 . With this in mind one can proceed to the path integral 
formulation, which gives origin to an integral over the Grassmann variables of an exponential con- 
taining a quadratic form in the fermionic fields, from H 2 and the normal ordering of He, as well as 
a quartic expression from : He :. The quartic expression can be reduced to a quadratic form by a 
Hubbard-Stratonovich transformation [10], through the introduction of a suitable auxiliary bosonic 
field (in our case a real field), and now the Gaussian integral over the fermionic variables can be 
explicitly performed, leaving an integral over the bosonic field only [11]. The problem, however, is 
to obtain an integral that can be interpreted as an integration over a well defined probabilistic mea- 
sure, which can be approximated by stochastic simulation techniques. In the following we show 
how this can be accomplished by taking advantage of the symmetries of the system 
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We start by rewriting the expression for the charge as 

q x = a\ x a x ^-a x -\a\_ v (4.1) 
We now introduce hole creation and annihilation operators for the electrons with spin — 1 : 

bl = a x -i, b x = al_ x (4.2) 

so that the charge becomes 

q x = ala x -blb x . (4.3) 

Note that we dropped the spin indices since from now on a,a* and b,b^ will always refer to spin 
1 and —1, respectively. Finally we change the sign of the b,b^ operators on one of the sublattices. 
The crucial constraint is that all redefinitions of the operators respect the anticommutator algebra. 
Because H2 only couples sites on the two different sublattices, it takes the form 

H 2 = £ -K(ala y + blb y + h.c). (4.4) 

(x,y) 

We introduce fermionic coherent states 

(VA*,T]*| = (0|e-^^^ + Mx) | V A,7]) = e -^(^4 + T ?^l)|0), (4.5) 

where \y x ,y*, ri x ,ri* are anticommuting fermionic variables (elements of a Grassmann algebra). 
The path integral formulation is obtained by factoring 

e~P H = e- H5 e~ H5 ...e~ H5 (N, terms) (4.6) 

with 8 = fi/N t , and then inserting repeatedly among the factors the resolution of the identity ex- 
pressed in terms of an integral over the fermionic variables. The trace in Eq. 3.4 must also be ex- 
pressed in terms of a similar integral (see e.g. [9] for details). This leads to integrals over fermionic 
fields Yx,t, Wxf> > Wxt ( tne index t = 0,---,N t — 1 appears because of the multiple resolutions 
of the identity and can be thought of as an index labeling Euclidean time), which contain in the 
integrand expressions of the type 

(vhXM~ H8 Wx,t,r} x ,t). (4.7) 

The last ingredient is the identity 

( ¥ :, t XAF(4M x ,a x ,b x )\ ¥x ^ (4.8) 

which is true of any normal ordered function F of the operators a\,b\,a x ,b x . 

As we indicated above, the Hamiltonian is in normal order form provided one separates the 
local term e 2 V xx q x q x in the Coulomb Hamiltonian into two normal-ordered pieces 

e 2 V xx q x q x = e 2 V xx : q x q x : + e 2 V xx {a\a x + b\b x ) . (4.9) 

By reassigning the quadratic term in Eq. 4.9 to H2, the exponent — H 8 in Eq. 4.7 is normal ordered. 
The exponential exp(— H 8) is not normal-ordered, but it differs from its normal ordered form by 
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terms 0(8 2 ), which give a vanishing contribution to the integral in the limit of N t — > °°. Neglecting 
these terms we may replace the operator expression exp(— H 8) with an exponential involving the 
fermionic fields, as follows from Eq. 4.8. We thus obtain the following expression for the partition 
function 

Z= lim /n^Y'm^n»>^ 

(4.10) 

where Q x t = Y xt ¥x,t — VxtVx^ and we have used m (and n) as a shorthand for the indices x,t. M is 
a matrix whose components may be deduced from 



£ WlMm,nWn = £ £ ¥xA¥x, t +l ~ Yx,t) + 5 - K £ (W*x,,¥y,t + Wy,tWx,t) 8 

m,n t x ( x .y) 



(4.11) 



where Xfajv, must be identified with —y/ x ,o- 

We now perform a Hubbard-Stratonovich transformation, introducing c-number real variables 
(p x .t to recast the exponential with the quartic term in the form 

e -Lx,y, t e 2 Qx,tVx,yQy, t 8 = f J~J ^ ^- Z x ,y, t $x,t ( V~ 1 )jt^S/4 g - Z x ,t ^x.tQxjS ^ ^ n) 

J x,t 

where we have absorbed a constant measure factor in the definition of the integral over (j> x>t . Insert- 
ing the r.h.s. of Eq. 4. 12 into Eq. 4. 10 and introducing the diagonal matrix & x ,t,y,T = (4>x,t8) 8x, y 8 t:T 
produces the compact result 

Z = J d<l>drdWn*dne-W^ 8 / 4 -Y*( M + ie ®)V-n*(M-ie®)r] (4 13) 

where we have used matrix notation for all the sums and have dropped the limit notation. 
The Gaussian integration over the anticommuting variables can now be done to obtain 

Z = J d(j)e~ ( l )V ' l( l )5 / 4 det(M-ie<^)det(M + ie^). (4.14) 

Because of the identity, 

det(M - ie<t>) det(M + ie<t>) = det[(M + ie<t>)\M + ie®)] 

the measure is positive definite. The down spins are treated as antiparticles (holes) moving back- 
ward in time relative to the up spins, exactly canceling the phase for each separately. Correla- 
tors for the fermion operators are now obtained by integrating the appropriate matrix elements of 
(M± ie<t>)~ 1 with the measure given by Eq. 4.14. 

Equation 4.14 is the main result of our work. It establishes the partition function and expecta- 
tion values as integrals over real variables with a positive definite measure. This is a crucial step for 
the application of stochastic approximation methods. There remains the problem of sampling the 
field (j) x>t with a measure which contains the determinant of a large matrix. But, following what is 
done in lattice gauge theory, this challenge can be overcome through the application of the hybrid 
Monte Carlo (HMC) technique [6]. In a broad outline, in HMC one first replaces the determinants 
in Eq. 4. 14 with a Gaussian integral over complex pseudofermionic variables £ x /. 

<kn[(M + w>t>?(M + u&)) = J d! ;* d t; e -£*(M+ie®y-\M+iecS>)- l C (4 15) 
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(In this equation and in the following Eq. 4.16 we absorb an irrelevant, constant measure factor in 
the integrals.) One then introduces real "momentum variables" % xf conjugate to <j> x>t and inserts in 
Eq. 4.15 unity written as a Gaussian integral over 71. One finally arrives at 

Z = |^^CVC^ V " 1 * a / 4 -^( M+ ^) t " 1 ( M + K *)"^-* 2 / 2 . (4-16) 

The idea of HMC is to consider the simultaneous distribution of the variables , 71, £ and £* deter- 
mined by the measure in Eq. 4. 16. The phase space of these variables is explored by first extracting 
the 71, £ and £* according to their Gaussian measure, and then evolving the and % variables with 
fixed £, according to the evolution determined by the Hamiltonian 

(tt,0) = ^ + +C*(M+^) t - 1 (M + ^)- 1 C. 

Because of Liouville's theorem, the combined motion through phase space produces an ensemble of 
variables distributed according to the measure in Eq. 4.16 and, in particular, of fields distributed 
according to Eq. 4.14. 

Of course, the discussion above assumes that the Hamiltonian evolution of and % is exact, 
which will not be the case with a numerical evolution. The HMC algorithm addresses this short- 
coming 1) by approximating the evolution with a symplectic integrator which is reversible and 
preserves phase space, 2) by performing a Metropolis accept-reject step at the end of the evolution, 
based on the variation of the value of the Hamiltonian. 

5. Numerical tests 

We tested our method on the two-site system obtained by taking L = 1 , which can be solved 
exactly. We label the sites x = 0, 1. With K = 1/3, the Hamiltonian H = Hi+Hq is now 

H 2 = -(ajao + «o a i +b\bo + blbi) +n(a*.a x + blb x ) 

2e 2 

He = 2e 2 (atao — blbo)(a\ai—b\bi)-\ cfi x b\a x b x (5.1) 

ro 

where we have taken Vb,i = V^o =1/3 and a local interaction term Vb,o = V^i = 1 /r$. The radius ro 
sets the physical scale in lattice units for localization of the net charge at the carbon atom. It must be 
restricted to ro < 1 for stability of the vacuum. Also the normal ordering prescription for e 2 V xx q x q x 
in Eq. 4.9 adds a new contribution to H2 in the form of an /3 "chemical potential" pLa\ s o^ s a x y . It 
is well known [12] that a 73 chemical potential for any value of /I does not introduce a phase in the 
measure. To maintain the full SU (2) "flavor" symmetry of the tight-binding graphene Hamiltonian, 
we must set ji. = e 2 /ro to its proper value. For the two-site system, the spin generators of Eq. 3.3 
become 

J + =A = {-Ifalbl and J 3 = \^ x a x + b\b x \j2 - 1, 

allowing us to unambiguously classify the 16 states as 5 singlets, 4 doublets and one triplet, given 
in Table 1. 
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Table 1: Spin multiplets and energies for the 2 site Graphene model. 



We compared KMC results for expectation values of several products of fermionic operators 
with the corresponding exact values, finding satisfactory agreement. For example, the correlation 
function 

C a {t) = ((ao- ai )(t)(al-al)(0))/2 (5.2) 

is illustrated in Figure. 5, which shows HMC results converging to the exact correlators for both 
the free theory with e = and an interacting case with e = 0.5. 




1 1 1 1 1 1_ 

U 1 2 3 4 5 * 

f 



Figure 5: HMC results for C a (t) (Eq. 5.2) with e = and e = 0.5 compared to the exact correlators. Here 
H = e 2 /r , r = 1/2 and J3 = N,8 = 6.4. 

A stringent test is to demonstrate the convergence to exact SU (2) symmetry in the "time" 
continuum limit. To this end, consider a second correlation function, 

C h (t) = ((&£ + &{)(/) (&o + &i)(0))/2 , (5.3) 
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\ i i 5 4 s 6 



Figure 6: HMC results for C a (t) (Eq. 5.2) and 
Q,(f) (Eq. 5.3) with e = 0.5 compared to the ex- 
act correlator. As in Figure. 5, jX — e 2 /ro, tq = \/2 
and/3 =N,8 = 6.4. 




■ut ■ ■ ■ i , , I i ■ ■ ■ .t 

n.mt n.ns hmo has n.M «.is 



Figure 7: Linear extrapolation (with error band) of 
the HMC energies for doublet correlators C a and 
Q, as a function of the "time" lattice spacing 8 = 
6A/N, for N, = 32, 64, 128 and 256. The dotted 
horizontal line marks Eq = 1.266. 



related to C a (t) by an St/ (2) rotation. Figure. 6 illustrates that HMC results for both C a {t) and 
Cb(t) converge to the same continuum limit. 

We extract the energies of the doublet states at nonzero 5 = P/N t by fitting the correlator 
data in Figure 6 to single exponentials, C(t) « e~ ^ t for fit range 0.4 < t < 4. The results in 
Figure. 7 clearly show linear behavior E « Z?o + ci 8, converging to the exact continuum Eq = e 2 + 
\/4 + e 4 — 1 « 1.266. The continuum limit is consistent with restoration of the SU (2) symmetry 
of the Hamiltonian: a joint linear fit to both sets of energies gives limg^o^ = 1-262 ±0.004, with 
c\ = 1.25 ±0.07 and —0.98 ±0.04 for correlators C a and Q,, respectively. 



6. Discussion 

Before scaling up our simulations to large lattices, it is useful to reflect on the formalism 
presented here. In our formulation, the theory is reduced to a 2+ 1 dimensional lattice theory with no 
sign problem. The non-local potential when expressed as a convolution can be computed efficiently 
using a Fast Fourier transform. The only approximation to the tight binding Hamiltonian for the 
nearest neighbor coupling is the discretization of the Euclidean time, 8. So as we extrapolate to 
8 — > 0, we can in principle solve the exact tight binding model numerically to arbitrary precession. 
However there are, as in all lattice field theories, improved discretization and algorithms that may 
accelerate the convergence. The difference here is the "continuum limit" involves only the "time" 
discretization and there is no symmetry rotating this axis into the hexagonal plane. So we are 
exploring alternative discretization in "time" which may improve computational efficiency and 
accuracy for small 8. 

On a theoretical level, it is instructive to re-derive the discretization for finite lattice spacing 8 
starting with our exact time continuum (8 = 0) acton or its Lagrangian (S = / cffJz?) derived above, 

&{t) = YU^(dt + ieo 3 ^ x (t))Y,{t) + e 2 V xx YUt)¥x(t) - * £ v4(0Vv(0 + t^CO^^CO , 

M 

(6.1) 
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where for simplicity we have combined the two spin components \\r, 77 into a spinor field: y/J = 
(~ty x ,f] x ) and y x = {^f x ,i] x ) T . Let us also introduce a a staggered mass parameter, +/ — m\j/ x Yx f° r 
sites on the A/B lattices respectively, which will be useful in numerical simulation in any case due 
to critical slowing down in the chiral limit. The result for the free theory (e = 0) on the nearest 
neighbor hexagonal lattice is to lift the zero mode symmetrically to give the continuum dispersion 
relation, 

E = ±i^co k 2 +m 2 . (6.2) 

At m = this is just the well know spectrum of the nearest neighbor graphene Hamiltonian H2 by 
going to Fourier space. Using the oblique reciprocal co-ordinates, k\ = 3k x /2 + y3k y /2,k2 = V3k y , 
on the Hexagonal lattice described in the caption of Figure. 6, (Ok is given by cot = | 1 + e^ 1 + e^ 2 | in 
spatial units set by taking K = 1. On a finite periodic lattice, the momenta are discrete (&, = 27inj/Lj 
for integer «,•), so to preserve the exact zero modes at the apex of the Dirac cones, k\ = — &2 = 
±27r/3, we must chose L; to be divisible by 3. 

Now let us re-introducing the temporal lattice spacing, 8, into our continuum Lagrangian 
(6.1). Our derivation above amounted to replacing the derivative by a forward difference, d t \j/ x — > 
{Wx.t+\ — Wx,t)/S, but there are other equally valid discretizations with the same continuum limit. 
Another attractive choice is to use this forward difference on the A lattice, d t \\f x — > (v^v,/+i — tyx,t) / S, 
but a backward difference on the B lattice, d t \\f x — > (Yx,t — Yx,t- which has an interesting lattice 
symmetry of time reversal times A/B lattice exchange. Comparing the free (e = 0) dispersion 
relation at finite 5, we have 

(2/8)e iE5 / 2 sin{E8/2) = ±i^co 2 +m 2 (6.3) 
for our original choice of all forward differences versus 



/ (O 2 -\- vn 2 

(2/5) sin(£5/2) = ±U . (6.4) 

y 1 — mo 

for the A/B alternating forward/ backward difference choice. Both dispersion relations of course 
yield the well know exact continuum dispersion relation (6.2) for 5—^0. However, the second 
one is accurate to 0(8 2 ) and has a symmetric spectrum in the complex plane characteristic of a 
relativistic fermion in the continuum. Very likely this form has some real advantages which we are 
actively investigating. 

Next let us return to the interacting case. Note in our Lagrangian (6.1), that the quadratic 
(kinetic) action for the electrons exhibits a 1-d space independent "gauge invariance" 

yr x (t) -> exp[-iea 3 d(t)]\j/x(t) , v4(0 -> </4(O ex P[^30(O] . <M0 -> <M0 +<^0(O , (6.5) 

much like the residual gauge invariance of the 4-d representation with gauge potential Aq. However 
is not the same object. The potential term in «Sf is not gauge invariant and thus it fixes this residual 
gauge transformation. It is natural to ask what happens if we preserve this gauge invariance for the 
kinetic term on the temporal lattice. We find an intriguing connection between the normal ordering 
term, e 2 V xx \l/ x Wx^ an d the gauge invariance of the kinetic term, which suggest another discretization 
using compact gauge variables on the links, 

J dtyl(t){d, + ieUt)^)Vx{t) -^£v4/(^' e5<T3 ^^+i " ¥x,t) , (6.6) 
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which like Wilson link variables in lattice field theory preserve gauge invariance of the kinetic term 
at the discrete times t for each temporal slice. Now let us expand in 8 as usual to see how this 
approaches the continuum Lagrangian (6.1 ): 

Y,VxAe ied ° 3 ^ xt Yx,t+i-¥x,t)^ J dt^dtYx + ie J dtylc^xWx- 8 J dtj(j) x 1 y}y x + - ■ ■ (6.7) 

The first and two second term on the right hand side are the required form to account for the 
Hubbard-Stratonovich transformation. Usually we would drop the third term as irrelevant in the 
continuum. However the Hubbard-Stratonovich field, <p x ,t, is not a dynamical field and on the lattice 
its fluctuations are controlled only by the quadratic potential in the action, 

1 f 8 

- J dt^V-^yit) -> £ - 0,Ay >.v/ - (6-8) 

This term implies that the variance for small 8 diverges on each time slice: (0 2 ) = 0(1/8). Con- 
sequently the third term is not suppressed and must be taken into account. Remarkably it can be 
shown that this term is exactly equivalent to adding the normal ordering term to leading order in 8. 
To show this we make a field redefinition: (j> x ~ <j> x — 8V xy <py\l/Jy y which when substituted into the 
lattice Lagrangian cancels the tp 2 term in Eq. 6.7 to leading order but to preserve the measure in the 
path integral also requires us to introduce the Jacobian of this change of variables. This Jacobian 
gives rise to our normal ordering term in the continuum limit, 



/ = J ®l, t det = J e Trl °g( l ~ 5V YV) ~ J ^ e S dtV ^(f)Vx(t) . 

(6.9) 

In summary an alternative to the normal ordering term on the lattice is to ignore it and introduce 
compact gauge links. 

An independent derivation that yields the same result is to reverse the argument in Sec. 4 
by first introducing the Hubbard-Stratonovich transformation in the Hamiltonian for each of the 
N t = jS/5factors in the discretized transfer matrix, 

e -pH = e -H8 e -H8 e -H8 ^ terms ). (6 10) 

Since the local charges commute {q x ,q y ] = the result of the Hubbard-Stratonovich transformation 
to leading order is simply the replacement, 



e 



-H8 



for each factor on time slice t. Now if we use the exact normal ordering identity, 

exp[-ie8(l) x ala 3 a x } -»• exp[y4e~ /e5a3< ^Vj , (6-12) 

given as Eq. A8 in Ref. [13] and expand to leading order, we again get the quadratic (seagull) term, 
— (e 2 /2)8 2 <p x 1 y%y x , in Eq. 6.7 . To formally take the continuum limit with this term of course 
we must again make a change of variables arriving to the unique continuum Lagrangian (6.1). 
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This alternate derivation makes it clear that the two lattice expressions (with — (e 2 /2)5 2 <^ l/zjl/^ or 
e 2 VxxWxWx) are simply different approximation to account for normal ordering. What is surprising 
is that compact link variables also provide this contribution and that in this form there is no explicit 
dependence on V xx need to define the effects of normal ordering. The dependence is implicit in 
the necessity to stabilize the quadratic form with a local repulsion at the Carbon ions to allow the 
Hubbard-Stratonovich transformation to be legitimate. 

Beyond the curiosity of this seagull term inducing the normal ordering correction, it opens up 
the attractive option of using compact links, as if were a gauge potential. One might suppose 
that this is a consequence of dimensional reduction from a real 4-d gauge theory but it is also more 
general since the result is exact for any static potential term V xy . The generality of the construction 
allows any static potential and any fixed geometry for the graphene sheet. It is also possible to in- 
clude fluctuation of the graphene sheet by modeling phonon interactions with additional dynamical 
spatial gauge links. At present we are investigating a range of discretization schemes numerically in 
comparison with solving the exact spectrum on lattices with L\ = L2 = 2 to see how each converges 
to the continuum limit. These test lattices have 8 sites with two spins per site so in total 2 16 = 65536 
states offer a rigorous and non-trivial verification of our HMC algorithm and software. 

7. Conclusion 

We have presented a Lagrangian formalism for the nearest neighbor tight binding theory for 
single layer graphene that allows one to do Hybrid Monte Carlos simulations by the traditional 
method of Euclidean lattice field theory. 

While the results we reported are for a small single hexagonal cell test systems, they demon- 
strate that HMC simulations of graphene directly on the hexagonal graphene lattice are possible and 
have the potential to produce valuable results. The dominant nearest neighbor hopping term has 
no sign problem, and we anticipate that a small next-to-nearest neighbor coupling k 1 /k ~ 0.03 can 
be accommodated by reweighting without a prohibitive cost. The crucial observation is the cancel- 
lation between the phase of the up spin and down spin determinant, when the latter are treated as 
holes moving backward in time. The only approximation is the discretization error introduce by a 
lattices spacing 8 for the "time" or temperature lattice. This opens up a method for effectively solv- 
ing this tight binding model numerically to arbitrary precision subject to taking the 1-d continuum 
limit. Relative to computational cost of lattice QCD simulations, clearly lattice graphene simula- 
tions are not only feasible but considerably more tractable with modern computing resources and 
algorithms. Still this is a young computational field and there is much work to be done to explore 
both algorithms and discretization scheme to optimize these simulations. 

We are currently pursuing simulations of larger systems, and beginning to explore the many 
possible generalizations such as distortions of the lattice, phonons, inclusion of magnetic fields, etc. 
We have a fully functioning code for our periodic hexagonal lattice with nearest neighbor hopping 
terms and a periodic Coulomb potential. We are putting the conjugate gradient inverter and the 
molecular dynamics evolution onto single GPUs. Since the NVIDIA Tesla class GPU have device 
memories in the range of 3 to 6 GigaBytes, this will allow a very substantial on card lattice volumes 
and a performance of 0(1 00) Gigaflops per GPU. The potential problems this open up to simulation 
are substantial. Graphene has a range of interesting properties. Our first goal is to repeat the 
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calculation in Ref. [4, 5] of the critical charge for the excitation gap (or chiral symmetry breaking) 
for single layer graphene due to the strong Coulomb potential. Having a reliable number for this 
critical parameter for the tight binding model is of considerable experimental interest. Beyond 
that there are a large range of effects due to changing boundary conditions and the introduction 
of phonons, which are easily introduced into the lattice Lagrangian formalism presented here. It 
is even feasible to simulate small graphene samples with lattices, identical in size and geometry 
to those being studied in experimental investigation. Here finite size effects are real and very 
interesting, both theoretically and technologically. 
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